stanEx4
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # not see parameters str1..., str2,... using regex as exc='[str1,str2,...]'
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}explanatory variable obs.x have noise
xN(x0.sx),yN(b0+b1*x0,s)
ex8-0.stan
normal regression without explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
normal regression with explanatory variable’s noise
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -88137.4
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 33 -37.0164 0.000318117 0.002379 1 1 69
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -37.02
## b0 7.54
## b1 1.74
## s 3.86
## m0[1] 33.62
## m0[2] 14.02
## m0[3] 37.27
## m0[4] 36.64
## m0[5] 37.17
## m0[6] 30.62
## m0[7] 31.04
## m0[8] 17.98
## m0[9] 8.31
## m0[10] 20.55
## m0[11] 9.64
## m0[12] 18.80
## m0[13] 14.36
## m0[14] 25.06
## m0[15] 20.36
## m0[16] 7.31
## m0[17] 18.85
## m0[18] 41.13
## m0[19] 32.61
## m0[20] 39.87
## y0[1] 26.14
## y0[2] 12.23
## y0[3] 38.67
## y0[4] 31.88
## y0[5] 37.67
## y0[6] 23.53
## y0[7] 33.60
## y0[8] 18.17
## y0[9] 13.15
## y0[10] 21.64
## y0[11] 10.44
## y0[12] 19.99
## y0[13] 7.96
## y0[14] 24.74
## y0[15] 19.88
## y0[16] 7.90
## y0[17] 16.70
## y0[18] 42.30
## y0[19] 36.45
## y0[20] 39.14
## m1[1] 7.31
## m1[2] 11.07
## m1[3] 14.82
## m1[4] 18.58
## m1[5] 22.34
## m1[6] 26.10
## m1[7] 29.85
## m1[8] 33.61
## m1[9] 37.37
## m1[10] 41.13
## y1[1] 7.16
## y1[2] 5.56
## y1[3] 15.33
## y1[4] 15.22
## y1[5] 24.53
## y1[6] 23.86
## y1[7] 27.92
## y1[8] 34.82
## y1[9] 34.94
## y1[10] 36.43
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -37.14 -36.86 1.19 1.00 -39.50 -35.84 1.01 625 753
## b0 7.60 7.55 1.66 1.61 4.83 10.37 1.01 348 421
## b1 1.74 1.74 0.15 0.15 1.49 1.97 1.01 505 691
## s 4.35 4.27 0.77 0.76 3.27 5.75 1.00 1498 1310
## m0[1] 33.62 33.65 1.24 1.18 31.57 35.58 1.00 2579 1537
## m0[2] 14.07 14.04 1.25 1.21 12.03 16.10 1.01 361 467
## m0[3] 37.26 37.28 1.46 1.41 34.82 39.56 1.00 1936 1210
## m0[4] 36.63 36.65 1.42 1.36 34.26 38.87 1.00 2058 1190
## m0[5] 37.16 37.18 1.46 1.40 34.74 39.45 1.00 1956 1213
## m0[6] 30.63 30.65 1.09 1.03 28.82 32.34 1.01 2697 1640
## m0[7] 31.05 31.07 1.11 1.05 29.21 32.79 1.00 2708 1674
## m0[8] 18.02 18.02 1.05 1.03 16.35 19.75 1.00 418 611
## m0[9] 8.38 8.33 1.61 1.56 5.69 11.02 1.01 347 420
## m0[10] 20.59 20.59 0.97 0.94 19.06 22.16 1.00 512 825
## m0[11] 9.70 9.66 1.52 1.46 7.19 12.19 1.01 348 430
## m0[12] 18.84 18.84 1.02 1.01 17.23 20.53 1.00 440 693
## m0[13] 14.41 14.39 1.23 1.20 12.40 16.42 1.01 363 486
## m0[14] 25.08 25.07 0.94 0.91 23.54 26.60 1.00 1012 1375
## m0[15] 20.39 20.39 0.98 0.95 18.86 21.97 1.00 502 811
## m0[16] 7.37 7.33 1.67 1.62 4.59 10.16 1.01 348 421
## m0[17] 18.88 18.89 1.02 1.01 17.28 20.58 1.00 441 703
## m0[18] 41.11 41.14 1.72 1.66 38.28 43.86 1.00 1448 1118
## m0[19] 32.61 32.64 1.19 1.12 30.66 34.48 1.00 2659 1639
## m0[20] 39.85 39.86 1.64 1.58 37.17 42.44 1.00 1576 1163
## y0[1] 33.57 33.63 4.58 4.56 25.90 40.92 1.00 1991 2045
## y0[2] 13.93 14.01 4.51 4.40 6.51 21.32 1.00 1411 1701
## y0[3] 37.22 37.30 4.60 4.25 29.67 44.45 1.00 2102 2089
## y0[4] 36.44 36.44 4.66 4.46 28.89 44.42 1.00 1844 1819
## y0[5] 37.00 37.12 4.73 4.52 29.16 44.30 1.00 2127 1933
## y0[6] 30.67 30.78 4.64 4.49 23.25 38.15 1.00 2137 2013
## y0[7] 31.06 30.93 4.54 4.30 23.53 38.60 1.00 2015 1855
## y0[8] 18.29 18.35 4.65 4.41 10.55 25.80 1.00 1725 1811
## y0[9] 8.29 8.31 4.85 4.77 0.32 16.00 1.00 1534 1786
## y0[10] 20.69 20.64 4.50 4.24 13.38 28.20 1.00 1716 1758
## y0[11] 9.84 9.81 4.64 4.62 2.56 17.41 1.00 1680 1859
## y0[12] 18.90 18.93 4.48 4.42 11.61 26.23 1.00 1634 1970
## y0[13] 14.57 14.46 4.58 4.25 7.09 22.20 1.00 1526 1529
## y0[14] 25.15 25.09 4.59 4.43 17.83 32.57 1.00 1854 1858
## y0[15] 20.37 20.37 4.44 4.18 13.21 27.69 1.00 1975 2035
## y0[16] 7.41 7.35 4.74 4.64 -0.40 15.04 1.00 1297 1487
## y0[17] 18.70 18.92 4.56 4.24 11.14 26.28 1.00 1719 1839
## y0[18] 41.05 41.00 4.77 4.52 33.11 48.95 1.00 1914 1974
## y0[19] 32.61 32.60 4.54 4.45 25.50 40.04 1.00 2097 1834
## y0[20] 39.77 39.90 4.71 4.44 31.53 47.29 1.00 2042 1971
## m1[1] 7.37 7.33 1.67 1.62 4.59 10.16 1.01 348 421
## m1[2] 11.12 11.09 1.42 1.38 8.79 13.44 1.01 350 426
## m1[3] 14.87 14.86 1.20 1.17 12.93 16.85 1.01 367 486
## m1[4] 18.62 18.62 1.03 1.01 17.00 20.32 1.00 434 693
## m1[5] 22.37 22.37 0.94 0.93 20.85 23.87 1.00 638 1006
## m1[6] 26.11 26.11 0.95 0.92 24.55 27.65 1.00 1253 1443
## m1[7] 29.86 29.87 1.06 1.00 28.11 31.54 1.01 2632 1668
## m1[8] 33.61 33.64 1.24 1.18 31.56 35.57 1.00 2580 1537
## m1[9] 37.36 37.38 1.47 1.41 34.92 39.67 1.00 1917 1210
## m1[10] 41.11 41.14 1.72 1.66 38.28 43.86 1.00 1448 1118
## y1[1] 7.28 7.21 4.57 4.44 -0.07 14.89 1.00 1302 1651
## y1[2] 11.17 11.13 4.68 4.61 3.88 18.91 1.00 1304 1698
## y1[3] 15.06 15.00 4.53 4.15 7.55 22.52 1.00 1499 1554
## y1[4] 18.71 18.79 4.54 4.18 11.22 26.27 1.00 1770 1912
## y1[5] 22.41 22.40 4.51 4.43 15.17 29.76 1.00 1756 1834
## y1[6] 26.17 26.06 4.61 4.52 18.61 33.80 1.00 1895 1788
## y1[7] 29.66 29.67 4.53 4.25 22.13 37.08 1.00 1959 1847
## y1[8] 33.69 33.68 4.50 4.10 26.42 41.18 1.00 2094 1727
## y1[9] 37.40 37.55 4.78 4.32 29.46 45.42 1.00 1881 2000
## y1[10] 41.19 41.18 4.74 4.51 33.39 48.82 1.00 1882 1892
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -44.46 -45.61 10.56 9.31 -59.89 -25.46 1.05 97 147
## b0 7.61 7.63 2.02 1.94 4.28 10.73 1.02 444 454
## b1 1.72 1.73 0.17 0.16 1.44 2.00 1.02 435 349
## s 2.95 2.98 1.39 1.57 0.74 5.15 1.02 125 85
## sx 1.77 1.79 0.89 0.88 0.36 3.22 1.02 132 116
## x0[1] 13.46 13.44 1.46 1.67 11.20 15.66 1.00 429 2495
## x0[2] 3.78 3.79 1.18 0.98 1.82 5.71 1.01 967 1025
## x0[3] 16.37 16.48 1.24 1.16 14.30 18.27 1.00 885 1793
## x0[4] 16.26 16.31 1.17 1.05 14.31 18.09 1.01 996 1930
## x0[5] 18.74 18.54 1.72 1.81 16.42 21.85 1.02 199 135
## x0[6] 13.27 13.27 1.14 0.99 11.36 15.02 1.01 747 1579
## x0[7] 15.69 15.57 1.89 2.25 13.07 18.90 1.02 176 167
## x0[8] 6.14 6.16 1.20 0.99 4.16 7.99 1.01 1400 1214
## x0[9] -1.40 -1.12 1.79 1.85 -4.73 0.91 1.02 214 172
## x0[10] 6.55 6.64 1.27 1.23 4.41 8.49 1.01 567 1543
## x0[11] 0.75 0.90 1.25 1.07 -1.43 2.61 1.01 606 1068
## x0[12] 7.53 7.46 1.35 1.42 5.48 9.65 1.01 275 668
## x0[13] 5.03 4.98 1.36 1.41 2.97 7.30 1.01 323 446
## x0[14] 8.73 8.78 1.40 1.48 6.45 10.82 1.01 416 1336
## x0[15] 9.11 9.05 1.65 1.91 6.81 11.70 1.01 221 763
## x0[16] -0.16 -0.10 1.38 1.09 -2.42 1.90 1.02 593 886
## x0[17] 7.23 7.18 1.23 1.20 5.35 9.21 1.01 480 644
## x0[18] 20.66 20.44 1.63 1.61 18.42 23.75 1.02 243 139
## x0[19] 13.54 13.59 1.25 1.22 11.52 15.43 1.00 694 2001
## x0[20] 17.38 17.43 1.35 1.46 15.21 19.37 1.00 564 1736
## m0[1] 30.83 30.66 2.62 2.94 27.04 35.20 1.01 299 1955
## m0[2] 14.19 14.13 2.00 1.68 11.02 17.55 1.01 2631 1856
## m0[3] 35.80 35.67 2.20 2.13 32.34 39.50 1.01 607 1679
## m0[4] 35.61 35.42 2.07 1.89 32.31 39.12 1.00 982 2456
## m0[5] 39.83 40.01 2.70 3.06 35.35 43.65 1.00 348 2469
## m0[6] 30.47 30.54 1.88 1.61 27.26 33.49 1.01 3531 1951
## m0[7] 34.59 34.66 2.99 3.70 29.82 38.98 1.01 246 2052
## m0[8] 18.23 18.23 1.94 1.66 15.09 21.34 1.01 2806 2005
## m0[9] 5.33 5.16 2.82 3.19 1.37 10.08 1.01 307 1717
## m0[10] 18.96 18.85 2.18 2.21 15.62 22.66 1.00 586 2779
## m0[11] 8.99 8.88 2.03 1.83 5.86 12.38 1.01 1677 1477
## m0[12] 20.61 20.72 2.23 2.34 16.94 24.05 1.01 351 2395
## m0[13] 16.33 16.52 2.27 2.39 12.54 19.77 1.01 389 1687
## m0[14] 22.71 22.59 2.38 2.61 19.12 26.52 1.01 354 2123
## m0[15] 23.31 23.35 2.73 3.16 19.00 27.22 1.01 265 1483
## m0[16] 7.45 7.50 2.10 1.84 3.96 10.77 1.01 2170 2309
## m0[17] 20.10 20.16 2.02 1.96 16.81 23.28 1.00 617 2276
## m0[18] 43.14 43.37 2.50 2.55 38.81 46.74 1.00 566 2375
## m0[19] 30.95 30.86 2.21 2.23 27.59 34.71 1.01 481 2026
## m0[20] 37.53 37.29 2.45 2.48 33.94 41.70 1.00 373 1977
## y0[1] 30.85 30.21 4.12 3.65 25.12 38.32 1.00 759 2570
## y0[2] 14.18 14.15 3.88 3.17 7.71 20.76 1.00 3469 3390
## y0[3] 35.81 35.41 4.01 3.29 29.70 42.66 1.01 2016 2266
## y0[4] 35.53 35.25 3.85 3.08 29.64 42.27 1.00 2742 2643
## y0[5] 39.83 40.52 4.21 3.68 32.17 45.71 1.00 822 2378
## y0[6] 30.51 30.52 3.77 3.07 24.15 36.92 1.01 3544 2486
## y0[7] 34.58 35.19 4.37 4.12 26.62 40.49 1.00 554 2235
## y0[8] 18.25 18.30 3.74 3.02 12.01 24.48 1.01 3994 2644
## y0[9] 5.30 4.60 4.32 3.90 -0.60 13.21 1.00 690 2048
## y0[10] 19.01 18.54 3.95 3.35 13.10 26.06 1.00 1806 2724
## y0[11] 8.97 8.76 3.81 3.11 2.84 15.49 1.01 3532 3162
## y0[12] 20.61 21.08 3.93 3.31 13.48 26.38 1.00 1532 3378
## y0[13] 16.26 16.79 4.07 3.45 8.86 22.17 1.00 1146 2783
## y0[14] 22.77 22.19 4.08 3.55 16.97 30.25 1.00 974 2469
## y0[15] 23.36 23.98 4.21 3.70 15.50 29.16 1.00 674 1684
## y0[16] 7.49 7.49 3.90 3.16 0.99 14.09 1.01 3634 2555
## y0[17] 20.13 20.40 3.85 3.26 13.66 26.26 1.00 2769 2631
## y0[18] 43.17 43.75 4.15 3.56 35.68 49.19 1.00 1300 2669
## y0[19] 30.95 30.46 3.94 3.23 25.29 38.08 1.00 1681 2195
## y0[20] 37.47 36.89 4.02 3.40 31.56 44.75 1.00 1154 1452
## m1[1] 7.39 7.40 2.04 1.96 4.03 10.53 1.02 442 454
## m1[2] 11.11 11.14 1.74 1.66 8.18 13.81 1.02 460 234
## m1[3] 14.83 14.87 1.46 1.33 12.35 17.17 1.02 490 218
## m1[4] 18.55 18.58 1.24 1.13 16.45 20.54 1.01 553 213
## m1[5] 22.27 22.27 1.09 1.01 20.43 24.06 1.01 620 215
## m1[6] 25.99 25.97 1.07 1.04 24.18 27.76 1.00 630 272
## m1[7] 29.71 29.69 1.16 1.13 27.84 31.68 1.01 638 611
## m1[8] 33.43 33.40 1.35 1.29 31.23 35.69 1.01 613 1146
## m1[9] 37.15 37.16 1.61 1.54 34.38 39.81 1.01 580 261
## m1[10] 40.87 40.88 1.90 1.81 37.59 44.00 1.01 552 234
## x1[1] -0.15 -0.14 1.95 1.43 -3.38 3.01 1.00 3829 3325
## x1[2] 2.01 2.01 2.03 1.49 -1.20 5.32 1.01 4205 2907
## x1[3] 4.15 4.15 1.94 1.50 0.95 7.31 1.01 3873 2222
## x1[4] 6.36 6.34 2.02 1.52 3.14 9.60 1.01 4070 1642
## x1[5] 8.44 8.45 1.91 1.45 5.35 11.58 1.01 4006 2326
## x1[6] 10.61 10.62 1.98 1.47 7.34 13.87 1.01 4055 2847
## x1[7] 12.75 12.80 1.98 1.48 9.44 15.95 1.01 4242 2669
## x1[8] 14.99 15.01 1.94 1.41 11.85 18.13 1.01 4156 2757
## x1[9] 17.12 17.14 1.97 1.49 13.91 20.32 1.01 3811 2789
## x1[10] 19.32 19.28 1.99 1.50 16.16 22.61 1.01 4114 2775
## y1[1] 7.42 7.43 3.81 3.40 1.35 13.55 1.00 1313 2054
## y1[2] 11.07 11.14 3.68 3.14 5.15 17.03 1.01 1569 2492
## y1[3] 14.82 14.80 3.48 2.95 9.09 20.67 1.00 1654 1962
## y1[4] 18.58 18.54 3.47 2.88 12.96 24.26 1.01 2646 2603
## y1[5] 22.34 22.23 3.48 2.84 16.79 28.36 1.00 2636 2514
## y1[6] 25.97 25.95 3.33 2.75 20.48 31.60 1.01 3566 3052
## y1[7] 29.67 29.56 3.42 2.82 24.18 35.30 1.00 3099 2663
## y1[8] 33.42 33.35 3.51 2.91 27.81 39.27 1.00 3098 2278
## y1[9] 37.16 37.10 3.56 3.07 31.42 43.01 1.00 2170 2859
## y1[10] 40.82 40.75 3.76 3.32 34.81 47.20 1.01 1710 3208
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
outlier
ex10.stan
objective variable have outlier, y~cauchy(b0+b1*x,s)
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -817.234
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -33.1259 0.000203143 0.00056136 1 1 28
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -33.13
## b0 6.30
## b1 1.90
## s 3.18
## m0[1] 12.01
## m0[2] 19.76
## m0[3] 14.92
## m0[4] 10.58
## m0[5] 6.38
## m0[6] 8.18
## m0[7] 16.78
## m0[8] 13.25
## m0[9] 12.50
## m0[10] 6.91
## m0[11] 8.62
## m0[12] 22.06
## m0[13] 20.27
## m0[14] 11.03
## m0[15] 15.72
## m0[16] 9.03
## m0[17] 22.80
## m0[18] 12.08
## m0[19] 17.88
## m0[20] 15.65
## y0[1] 10.05
## y0[2] 20.11
## y0[3] 10.98
## y0[4] 8.94
## y0[5] 5.57
## y0[6] 9.06
## y0[7] 15.88
## y0[8] 13.73
## y0[9] 8.03
## y0[10] 0.33
## y0[11] 7.52
## y0[12] 20.30
## y0[13] 20.03
## y0[14] 10.39
## y0[15] 16.68
## y0[16] 17.01
## y0[17] 19.68
## y0[18] 12.49
## y0[19] 13.15
## y0[20] 13.93
## m1[1] 6.38
## m1[2] 8.21
## m1[3] 10.03
## m1[4] 11.86
## m1[5] 13.68
## m1[6] 15.51
## m1[7] 17.33
## m1[8] 19.15
## m1[9] 20.98
## m1[10] 22.80
## y1[1] 8.43
## y1[2] 5.43
## y1[3] 7.78
## y1[4] 9.01
## y1[5] 13.08
## y1[6] 9.33
## y1[7] 17.45
## y1[8] 16.43
## y1[9] 24.82
## y1[10] 27.25
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.55 -33.26 1.27 1.10 -36.15 -32.15 1.00 489 768
## b0 6.24 6.33 1.52 1.35 3.62 8.76 1.01 280 298
## b1 1.92 1.91 0.32 0.30 1.40 2.45 1.00 385 486
## s 3.60 3.51 0.65 0.60 2.71 4.77 1.00 1410 1402
## m0[1] 11.98 11.98 0.88 0.83 10.50 13.43 1.00 368 569
## m0[2] 19.79 19.73 1.26 1.19 17.68 21.88 1.00 1292 1424
## m0[3] 14.91 14.90 0.83 0.78 13.59 16.31 1.00 962 984
## m0[4] 10.55 10.58 0.99 0.95 8.90 12.19 1.01 308 440
## m0[5] 6.32 6.42 1.51 1.34 3.73 8.83 1.01 280 303
## m0[6] 8.13 8.20 1.27 1.17 6.02 10.23 1.01 281 326
## m0[7] 16.78 16.76 0.93 0.87 15.23 18.33 1.00 1757 1280
## m0[8] 13.23 13.24 0.82 0.79 11.93 14.61 1.00 498 789
## m0[9] 12.48 12.49 0.85 0.81 11.06 13.88 1.00 408 654
## m0[10] 6.85 6.95 1.44 1.29 4.40 9.26 1.01 280 307
## m0[11] 8.57 8.64 1.21 1.14 6.58 10.59 1.01 283 326
## m0[12] 22.10 22.04 1.57 1.51 19.48 24.72 1.00 959 1168
## m0[13] 20.30 20.25 1.32 1.26 18.10 22.52 1.00 1194 1382
## m0[14] 11.00 11.01 0.95 0.91 9.38 12.59 1.00 322 478
## m0[15] 15.72 15.70 0.86 0.78 14.33 17.16 1.00 1362 1150
## m0[16] 8.99 9.05 1.16 1.10 7.07 10.94 1.01 286 353
## m0[17] 22.85 22.80 1.68 1.63 20.05 25.68 1.00 892 1135
## m0[18] 12.05 12.06 0.87 0.84 10.57 13.50 1.00 373 585
## m0[19] 17.90 17.86 1.04 0.96 16.15 19.65 1.00 1714 1165
## m0[20] 15.65 15.63 0.86 0.78 14.27 17.08 1.00 1330 1150
## y0[1] 12.02 12.03 3.78 3.61 5.90 18.18 1.00 1736 1998
## y0[2] 19.66 19.74 3.83 3.68 13.40 25.92 1.00 1673 1654
## y0[3] 14.95 14.96 3.69 3.48 8.75 20.85 1.00 1772 1853
## y0[4] 10.40 10.35 3.94 3.85 3.84 16.95 1.00 1733 1691
## y0[5] 6.38 6.35 4.04 4.07 -0.19 12.98 1.00 1028 1643
## y0[6] 8.10 8.06 3.80 3.60 1.86 14.48 1.00 1177 1844
## y0[7] 16.78 16.88 3.81 3.66 10.43 22.67 1.00 1958 1860
## y0[8] 13.25 13.23 3.80 3.60 6.75 19.33 1.00 1950 1896
## y0[9] 12.54 12.40 3.72 3.52 6.56 18.93 1.00 1762 1974
## y0[10] 6.90 6.90 3.86 3.73 0.49 13.18 1.00 957 1642
## y0[11] 8.45 8.41 3.70 3.52 2.45 14.40 1.00 1531 1807
## y0[12] 21.96 22.03 3.99 3.90 15.42 28.61 1.00 1826 1740
## y0[13] 20.44 20.39 3.84 3.80 14.16 26.78 1.00 1874 2060
## y0[14] 10.83 10.82 3.70 3.66 4.86 16.75 1.00 1342 1884
## y0[15] 15.74 15.70 3.84 3.77 9.66 21.90 1.00 2070 1974
## y0[16] 8.91 8.89 3.83 3.67 2.72 15.33 1.00 1387 1815
## y0[17] 22.80 22.74 4.06 3.89 15.97 29.51 1.00 1834 1942
## y0[18] 12.14 12.12 3.78 3.62 5.81 18.35 1.00 1339 1695
## y0[19] 17.99 18.01 3.75 3.36 11.65 24.24 1.00 2035 1800
## y0[20] 15.53 15.40 3.85 3.74 9.49 21.99 1.00 1998 1920
## m1[1] 6.32 6.42 1.51 1.34 3.73 8.83 1.01 280 303
## m1[2] 8.16 8.23 1.26 1.17 6.06 10.26 1.01 281 324
## m1[3] 9.99 10.03 1.05 0.99 8.28 11.75 1.01 297 423
## m1[4] 11.83 11.83 0.89 0.85 10.34 13.31 1.00 359 562
## m1[5] 13.67 13.68 0.81 0.76 12.35 15.03 1.00 576 841
## m1[6] 15.50 15.48 0.85 0.77 14.13 16.92 1.00 1250 1084
## m1[7] 17.34 17.31 0.98 0.91 15.70 19.00 1.00 1796 1213
## m1[8] 19.18 19.13 1.18 1.11 17.20 21.19 1.00 1425 1306
## m1[9] 21.01 20.95 1.42 1.34 18.66 23.40 1.00 1087 1354
## m1[10] 22.85 22.80 1.68 1.63 20.05 25.68 1.00 892 1135
## y1[1] 6.31 6.24 3.92 3.71 -0.14 12.75 1.00 1020 1533
## y1[2] 8.12 8.25 3.83 3.77 1.74 14.24 1.00 1805 1803
## y1[3] 10.08 10.06 3.85 3.75 3.70 16.22 1.00 1451 1774
## y1[4] 11.88 11.88 3.71 3.49 5.78 17.77 1.00 1726 1851
## y1[5] 13.62 13.63 3.86 3.68 7.16 20.01 1.00 1789 1919
## y1[6] 15.53 15.59 3.71 3.66 9.49 21.62 1.00 1734 1969
## y1[7] 17.27 17.29 3.80 3.57 11.07 23.44 1.00 1919 1693
## y1[8] 19.16 19.23 3.83 3.64 12.91 25.29 1.00 2038 2089
## y1[9] 20.98 21.06 4.03 3.89 14.29 27.52 1.00 1827 1872
## y1[10] 22.94 22.92 4.02 4.04 16.48 29.41 1.00 1643 1810
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -101.724
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -14.165 0.000738009 0.00153556 0.9772 0.9772 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -14.16
## b0 5.15
## b1 2.04
## s 0.49
## m0[1] 11.28
## m0[2] 19.59
## m0[3] 14.40
## m0[4] 9.75
## m0[5] 5.25
## m0[6] 7.17
## m0[7] 16.39
## m0[8] 12.61
## m0[9] 11.81
## m0[10] 5.80
## m0[11] 7.64
## m0[12] 22.05
## m0[13] 20.14
## m0[14] 10.23
## m0[15] 15.25
## m0[16] 8.09
## m0[17] 22.85
## m0[18] 11.35
## m0[19] 17.58
## m0[20] 15.19
## y0[1] 11.57
## y0[2] 19.59
## y0[3] 14.87
## y0[4] 9.36
## y0[5] 11.07
## y0[6] 4.19
## y0[7] 14.81
## y0[8] 11.15
## y0[9] 43.30
## y0[10] 5.98
## y0[11] 7.72
## y0[12] 23.50
## y0[13] 20.67
## y0[14] 10.31
## y0[15] 15.04
## y0[16] 8.02
## y0[17] 20.20
## y0[18] 11.54
## y0[19] 17.53
## y0[20] 14.52
## m1[1] 5.25
## m1[2] 7.20
## m1[3] 9.16
## m1[4] 11.11
## m1[5] 13.07
## m1[6] 15.03
## m1[7] 16.98
## m1[8] 18.94
## m1[9] 20.90
## m1[10] 22.85
## y1[1] 10.20
## y1[2] 7.50
## y1[3] 7.32
## y1[4] 11.18
## y1[5] 12.88
## y1[6] 14.78
## y1[7] 14.36
## y1[8] 21.08
## y1[9] 18.71
## y1[10] 23.10
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.63 -16.32 1.40 1.24 -19.31 -15.04 1.00 438 597
## b0 5.21 5.19 0.37 0.30 4.65 5.84 1.01 566 482
## b1 2.03 2.04 0.07 0.06 1.91 2.14 1.01 789 743
## s 0.74 0.71 0.26 0.24 0.40 1.21 1.01 488 503
## m0[1] 11.31 11.29 0.24 0.21 10.92 11.71 1.01 711 801
## m0[2] 19.58 19.58 0.31 0.28 19.07 20.08 1.01 2054 1272
## m0[3] 14.41 14.41 0.23 0.22 14.04 14.77 1.00 1292 1146
## m0[4] 9.79 9.77 0.26 0.23 9.37 10.22 1.01 615 648
## m0[5] 5.30 5.28 0.37 0.30 4.74 5.93 1.01 565 482
## m0[6] 7.21 7.19 0.32 0.27 6.73 7.76 1.01 562 539
## m0[7] 16.40 16.40 0.25 0.23 16.00 16.80 1.01 1981 1038
## m0[8] 12.63 12.63 0.23 0.21 12.26 13.01 1.01 868 919
## m0[9] 11.83 11.82 0.23 0.21 11.45 12.23 1.01 762 878
## m0[10] 5.86 5.83 0.35 0.30 5.33 6.46 1.01 562 481
## m0[11] 7.69 7.67 0.31 0.26 7.22 8.21 1.01 566 544
## m0[12] 22.03 22.03 0.38 0.33 21.42 22.64 1.00 1746 1245
## m0[13] 20.13 20.13 0.32 0.29 19.60 20.66 1.01 1982 1312
## m0[14] 10.26 10.24 0.25 0.23 9.86 10.68 1.01 638 706
## m0[15] 15.26 15.27 0.23 0.22 14.89 15.64 1.00 1580 1000
## m0[16] 8.13 8.11 0.30 0.25 7.68 8.63 1.01 572 600
## m0[17] 22.83 22.83 0.40 0.35 22.18 23.47 1.00 1659 1311
## m0[18] 11.38 11.37 0.24 0.21 11.00 11.78 1.01 717 801
## m0[19] 17.58 17.58 0.27 0.24 17.14 18.01 1.01 2146 1075
## m0[20] 15.20 15.20 0.23 0.22 14.82 15.58 1.00 1558 1004
## y0[1] 11.88 11.29 22.52 1.09 6.93 16.68 1.00 1826 2015
## y0[2] 19.65 19.57 12.01 1.08 14.51 24.55 1.00 1921 1430
## y0[3] 13.51 14.41 36.49 1.06 9.60 18.42 1.00 1829 1988
## y0[4] 10.31 9.77 31.16 1.14 4.66 14.90 1.00 1888 1878
## y0[5] 4.62 5.30 20.37 1.12 -0.09 10.11 1.00 1869 1927
## y0[6] 6.17 7.17 44.66 1.08 2.32 11.65 1.00 2025 2055
## y0[7] 18.84 16.36 87.97 1.07 11.21 20.85 1.00 1834 1832
## y0[8] 12.80 12.65 24.53 1.03 7.65 17.02 1.00 1851 2011
## y0[9] 12.20 11.83 25.40 1.10 6.00 16.77 1.00 2077 1971
## y0[10] 7.18 5.83 44.20 1.20 0.51 11.61 1.00 1995 2042
## y0[11] 6.17 7.66 58.92 1.11 2.68 11.82 1.00 1911 1882
## y0[12] 21.37 22.02 18.48 1.13 17.13 26.02 1.00 1908 1847
## y0[13] 17.74 20.12 528.57 1.15 15.64 25.02 1.00 1822 1840
## y0[14] 8.35 10.24 93.64 1.05 6.23 15.08 1.00 1912 1885
## y0[15] 15.57 15.28 27.93 1.11 10.66 20.60 1.00 1993 2057
## y0[16] 8.06 8.11 18.39 1.11 3.69 12.47 1.00 2039 1855
## y0[17] 23.06 22.88 24.58 1.15 18.58 27.78 1.00 1996 2058
## y0[18] 12.30 11.37 32.48 1.10 6.81 15.95 1.00 1953 1786
## y0[19] 18.16 17.61 38.75 1.16 13.13 22.46 1.00 2076 2010
## y0[20] 6.71 15.19 324.96 1.14 10.81 19.74 1.00 2114 1972
## m1[1] 5.30 5.28 0.37 0.30 4.74 5.93 1.01 565 482
## m1[2] 7.25 7.23 0.32 0.27 6.77 7.79 1.01 563 558
## m1[3] 9.20 9.18 0.27 0.24 8.77 9.65 1.01 593 611
## m1[4] 11.14 11.13 0.24 0.22 10.76 11.55 1.01 699 762
## m1[5] 13.09 13.09 0.23 0.22 12.73 13.46 1.00 947 1016
## m1[6] 15.04 15.04 0.23 0.22 14.67 15.42 1.00 1499 1047
## m1[7] 16.99 16.99 0.26 0.24 16.57 17.41 1.00 2092 1000
## m1[8] 18.93 18.94 0.30 0.27 18.45 19.42 1.01 2124 1292
## m1[9] 20.88 20.89 0.34 0.31 20.32 21.44 1.00 1888 1289
## m1[10] 22.83 22.83 0.40 0.35 22.18 23.47 1.00 1659 1311
## y1[1] 4.85 5.28 33.30 1.20 0.24 10.33 1.00 1702 1910
## y1[2] 7.61 7.23 24.64 1.17 2.56 13.01 1.00 1965 1857
## y1[3] 9.27 9.17 30.60 1.02 5.12 13.33 1.00 1472 1594
## y1[4] 7.73 11.11 109.98 1.12 6.14 15.54 1.00 1901 1932
## y1[5] 16.21 13.07 134.29 1.03 8.68 17.60 1.00 1853 2008
## y1[6] 16.15 15.02 26.11 1.06 10.87 19.56 1.00 1947 1854
## y1[7] 17.52 17.01 24.32 1.10 12.61 22.50 1.00 1678 1738
## y1[8] 18.90 18.88 18.82 1.14 13.96 23.59 1.00 2130 1940
## y1[9] 19.14 20.88 126.44 1.12 16.32 26.22 1.00 2237 2039
## y1[10] 24.61 22.86 45.36 1.24 18.23 27.63 1.00 1964 1692
censored data
objective variable has NA
(x,y) i=1-N
(x0,y0) i=1-N0
x1 i=1-N1, y1=NA
(x,y)~N((mx,my),(sx2,sy2,sxy))
(x0,y0)~N((mx,my),(sx2,sy2,sxy))
x1~N(mx,sx2)ex11-0.stan
data{
int N0;
array[N0] vector[2] xy;
int N1;
vector[N1] x1;
}
parameters{
vector[2] m;
cov_matrix[2] s;
}
model{
target+=multi_normal_lpdf(xy | m, s);
x1~normal(m[1],s[1,1]^.5);
}
generated quantities{
vector[2] xy1;
xy1=multi_normal_rng(m,s);
real cr;
cr=s[1,2]/(s[1,1]*s[2,2])^.5;
}n=30
x=runif(n,0,9)
y=rnorm(n,10+3*x,4)
cor(x,y)## [1] 0.926qplot(x,y)
L=4
n0=sum(x>L)
x0=x[x>L]
y0=y[x>L]
x1=x[x<=L]
n1=sum(x<=L)
cor(x0,y0)## [1] 0.837qplot(x0,y0)
mdl=cmdstan_model('./ex11-0.stan')
data=list(N0=n0,xy=matrix(c(x0,y0),ncol=2),N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -8860.57
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 50 -100.084 0.00102871 0.00225021 1 1 77
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -100.08
## m[1] 4.25
## m[2] 19.37
## s[1,1] 6.81
## s[2,1] 24.16
## s[1,2] 24.16
## s[2,2] 96.37
## xy1[1] 3.39
## xy1[2] 15.39
## cr 0.94mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.4 seconds.
## Chain 2 finished in 0.4 seconds.
## Chain 3 finished in 0.4 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.4 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.4 seconds.
## Total execution time: 0.6 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -96.17 -95.81 1.75 1.63 -99.56 -94.00 1.01 507 815
## m[1] 4.27 4.28 0.55 0.55 3.37 5.21 1.01 391 388
## m[2] 19.48 19.62 2.76 2.79 15.04 23.76 1.02 201 115
## s[1,1] 8.79 8.35 2.69 2.33 5.34 13.89 1.00 814 943
## s[2,1] 30.80 29.26 11.24 9.99 15.39 51.53 1.02 265 433
## s[1,2] 30.80 29.26 11.24 9.99 15.39 51.53 1.02 265 433
## s[2,2] 129.97 118.95 57.31 51.11 56.10 235.03 1.02 210 397
## xy1[1] 4.28 4.30 2.89 2.81 -0.49 8.89 1.00 1915 1876
## xy1[2] 19.42 19.91 11.17 10.74 1.22 36.06 1.00 1497 833
## cr 0.92 0.93 0.06 0.04 0.81 0.97 1.02 476 523
xy=mcmc$draws('xy1')
cor(xy[,,1],xy[,,2])## [1] 0.889qplot(xy[,,1],xy[,,2])
objective variable is censored
y i=1-N, y~N(m,s)
actual ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11-1.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitionmdl=cmdstan_model('./ex8-0.stan')
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mle=mdl$optimize(data=data)## Initial log joint probability = -20850.3
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 25 -15.3216 0.0011708 0.000321661 1 1 55
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -15.32
## b0 17.34
## b1 1.34
## s 3.33
## m0[1] 25.91
## m0[2] 21.57
## m0[3] 22.06
## m0[4] 23.04
## m0[5] 21.61
## m0[6] 26.18
## m0[7] 24.32
## m0[8] 23.91
## m0[9] 24.67
## y0[1] 26.97
## y0[2] 19.61
## y0[3] 23.36
## y0[4] 22.76
## y0[5] 25.05
## y0[6] 26.71
## y0[7] 22.03
## y0[8] 24.08
## y0[9] 20.02
## m1[1] 18.89
## m1[2] 20.06
## m1[3] 21.23
## m1[4] 22.39
## m1[5] 23.56
## m1[6] 24.73
## m1[7] 25.89
## m1[8] 27.06
## m1[9] 28.22
## m1[10] 29.39
## y1[1] 16.50
## y1[2] 23.37
## y1[3] 15.09
## y1[4] 26.74
## y1[5] 25.46
## y1[6] 26.43
## y1[7] 28.54
## y1[8] 28.90
## y1[9] 30.68
## y1[10] 31.85mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.97 -15.55 1.62 1.26 -19.02 -14.31 1.02 356 362
## b0 18.15 18.17 6.75 6.22 7.99 28.76 1.04 281 213
## b1 1.19 1.19 1.39 1.22 -1.02 3.25 1.04 310 231
## s 4.73 4.33 1.82 1.27 2.84 7.83 1.00 843 836
## m0[1] 25.74 25.77 2.90 2.45 21.05 30.08 1.01 625 521
## m0[2] 21.90 21.85 2.75 2.39 17.70 26.19 1.02 296 314
## m0[3] 22.33 22.28 2.37 2.07 18.66 26.07 1.02 321 343
## m0[4] 23.20 23.14 1.82 1.52 20.33 26.13 1.01 512 707
## m0[5] 21.93 21.88 2.71 2.35 17.77 26.18 1.02 298 299
## m0[6] 25.98 26.03 3.14 2.63 20.93 30.66 1.01 571 510
## m0[7] 24.33 24.35 1.85 1.54 21.33 27.22 1.00 1249 981
## m0[8] 23.96 23.94 1.73 1.41 21.24 26.72 1.00 1274 972
## m0[9] 24.64 24.66 2.01 1.66 21.48 27.73 1.00 1100 850
## y0[1] 25.77 25.87 5.74 5.07 16.74 34.58 1.00 1474 1344
## y0[2] 21.84 21.65 5.49 4.93 13.50 30.89 1.00 762 1209
## y0[3] 22.40 22.21 5.80 4.67 13.94 31.49 1.00 1047 1130
## y0[4] 23.38 23.27 5.51 4.53 15.27 31.62 1.00 1646 1785
## y0[5] 21.93 21.85 5.80 4.86 13.08 31.22 1.01 879 1301
## y0[6] 25.76 25.66 6.05 5.04 16.37 35.26 1.00 1268 1215
## y0[7] 24.30 24.29 5.46 4.76 15.97 33.29 1.00 2003 1541
## y0[8] 24.14 24.04 5.42 4.68 15.86 32.94 1.00 1778 1535
## y0[9] 24.61 24.58 5.46 4.82 16.06 33.59 1.00 1859 1891
## m1[1] 19.52 19.54 5.21 4.73 11.83 27.88 1.04 278 213
## m1[2] 20.56 20.56 4.09 3.68 14.49 27.09 1.04 278 222
## m1[3] 21.59 21.52 3.03 2.66 17.00 26.30 1.03 288 260
## m1[4] 22.62 22.57 2.15 1.89 19.27 25.98 1.01 354 387
## m1[5] 23.66 23.61 1.71 1.38 20.94 26.42 1.00 851 931
## m1[6] 24.69 24.71 2.04 1.67 21.44 27.81 1.00 1072 849
## m1[7] 25.72 25.76 2.89 2.43 21.05 30.04 1.01 628 521
## m1[8] 26.76 26.83 3.93 3.29 20.42 32.59 1.01 473 453
## m1[9] 27.79 27.93 5.04 4.26 19.87 35.19 1.01 415 382
## m1[10] 28.82 28.94 6.19 5.24 19.01 38.01 1.02 385 346
## y1[1] 19.46 19.53 7.56 6.23 7.83 31.33 1.01 413 602
## y1[2] 20.51 20.37 6.43 5.52 9.92 30.74 1.01 558 747
## y1[3] 21.60 21.47 6.02 5.31 12.21 30.95 1.00 853 1221
## y1[4] 22.69 22.55 5.29 4.67 14.20 31.03 1.00 896 1095
## y1[5] 23.55 23.57 5.31 4.54 15.22 31.71 1.00 1828 1599
## y1[6] 24.73 24.62 5.57 4.66 15.94 33.97 1.00 1788 1572
## y1[7] 25.90 25.98 5.89 4.81 16.41 34.85 1.00 1442 1253
## y1[8] 26.75 26.74 6.47 5.42 16.71 37.25 1.00 950 927
## y1[9] 28.02 28.25 7.28 6.18 16.91 39.39 1.01 707 722
## y1[10] 28.68 28.76 7.93 6.74 15.87 41.41 1.01 541 474
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11-1.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -333.417
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 25 -21.8357 0.000501952 4.20973e-05 0.9955 0.9955 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -21.84
## b0 5.78
## b1 3.93
## s 4.51
## m0[1] 30.94
## m0[2] 18.21
## m0[3] 19.64
## m0[4] 22.53
## m0[5] 18.33
## m0[6] 31.75
## m0[7] 26.29
## m0[8] 25.06
## m0[9] 27.30
## y0[1] 29.13
## y0[2] 21.85
## y0[3] 19.29
## y0[4] 22.04
## y0[5] 19.63
## y0[6] 32.25
## y0[7] 28.27
## y0[8] 29.41
## y0[9] 24.75
## m1[1] 10.35
## m1[2] 13.77
## m1[3] 17.20
## m1[4] 20.62
## m1[5] 24.05
## m1[6] 27.47
## m1[7] 30.90
## m1[8] 34.32
## m1[9] 37.74
## m1[10] 41.17
## y1[1] 16.86
## y1[2] 10.31
## y1[3] 9.99
## y1[4] 27.31
## y1[5] 20.55
## y1[6] 29.08
## y1[7] 33.12
## y1[8] 38.61
## y1[9] 43.36
## y1[10] 39.91mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.08 -21.72 1.48 1.17 -24.85 -20.52 1.02 325 546
## b0 1.58 2.45 6.69 5.76 -11.81 10.30 1.04 105 94
## b1 4.88 4.64 1.42 1.20 3.06 7.50 1.04 100 104
## s 6.42 5.91 2.43 1.81 3.73 10.89 1.02 202 499
## m0[1] 32.81 32.25 3.44 2.66 28.65 38.72 1.02 204 286
## m0[2] 17.01 17.41 2.78 2.35 11.60 20.89 1.02 199 271
## m0[3] 18.79 19.06 2.46 2.10 14.25 22.35 1.02 283 371
## m0[4] 22.37 22.39 2.09 1.75 19.06 25.47 1.01 701 743
## m0[5] 17.15 17.53 2.75 2.32 11.83 21.00 1.02 203 293
## m0[6] 33.82 33.20 3.68 2.83 29.43 40.24 1.02 182 246
## m0[7] 27.04 26.77 2.32 1.80 24.07 30.69 1.01 1429 716
## m0[8] 25.52 25.34 2.16 1.74 22.52 28.85 1.01 1410 965
## m0[9] 28.29 27.98 2.51 1.95 25.18 32.36 1.01 738 619
## y0[1] 32.73 32.09 7.72 6.28 21.35 45.22 1.01 1541 1105
## y0[2] 17.28 17.59 7.34 6.48 5.18 28.32 1.01 1174 1388
## y0[3] 19.01 19.09 7.16 6.17 7.63 30.39 1.00 1685 1257
## y0[4] 22.56 22.69 7.17 6.32 10.85 34.13 1.00 1714 1165
## y0[5] 17.42 17.76 7.43 6.35 5.16 28.60 1.01 1102 1042
## y0[6] 33.86 33.35 7.62 6.52 22.74 46.92 1.01 605 827
## y0[7] 26.97 26.77 7.16 6.34 15.89 38.60 1.00 1764 1805
## y0[8] 25.56 25.48 7.18 6.09 14.37 36.87 1.00 1807 1545
## y0[9] 28.27 28.04 7.12 6.16 17.37 40.04 1.00 1825 1508
## m1[1] 7.25 8.01 5.14 4.40 -3.19 14.04 1.03 110 88
## m1[2] 11.50 12.14 4.04 3.49 3.64 16.84 1.03 125 158
## m1[3] 15.75 16.21 3.04 2.62 9.85 19.84 1.03 171 248
## m1[4] 20.00 20.18 2.29 1.93 15.94 23.33 1.01 359 433
## m1[5] 24.25 24.15 2.08 1.70 21.19 27.31 1.01 1186 991
## m1[6] 28.50 28.17 2.55 1.98 25.36 32.65 1.01 659 564
## m1[7] 32.75 32.20 3.42 2.64 28.60 38.64 1.02 206 286
## m1[8] 37.00 36.22 4.47 3.48 31.63 44.99 1.03 146 199
## m1[9] 41.25 40.36 5.60 4.43 34.57 51.36 1.03 126 124
## m1[10] 45.50 44.33 6.77 5.43 37.34 57.69 1.03 116 122
## y1[1] 7.28 8.02 8.66 7.72 -7.80 19.34 1.01 453 388
## y1[2] 11.53 12.00 7.93 6.73 -2.11 22.87 1.00 803 956
## y1[3] 15.83 16.08 7.47 6.54 2.89 27.00 1.00 972 1512
## y1[4] 19.98 20.15 7.53 6.39 7.90 31.54 1.01 1708 1280
## y1[5] 24.18 24.20 7.33 6.58 12.53 35.73 1.00 1675 1698
## y1[6] 28.26 27.87 7.52 6.40 17.20 40.63 1.00 1935 1364
## y1[7] 32.75 32.16 7.67 6.60 21.23 45.45 1.01 1195 1361
## y1[8] 37.11 36.26 8.27 7.00 25.65 51.35 1.00 1141 1004
## y1[9] 41.11 40.17 8.99 7.32 28.91 57.03 1.01 387 644
## y1[10] 45.18 44.09 9.28 7.80 32.69 60.56 1.02 256 437
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
sensitivity/specificity analysis
ex14.stan
estimating sens and spec
data {
int N;
array[N] int x;
array[N] int y;
}
parameters {
real<lower=0,upper=1> p;
real<lower=0,upper=1> se;
real<lower=0,upper=1> sp;
}
model {
p~uniform(0,1);
se~uniform(0,1);
sp~uniform(0,1);
for (i in 1:N) {
y[i]~bernoulli(x[i]*se+(1-x[i])*(1-sp));
}
}
generated quantities {
real ppv;
real npv;
ppv=se*p/((se*p)+((1-p)*(1-sp)));
npv=(1-p)*sp/(((1-p)*sp)+(p*(1-se)));
}n=20
data=list(N=n,
x=sample(0:1,n,replace=T),
y=sample(0:1,n,replace=T))
mdl=cmdstan_model('./ex14.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -11.4166
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -11.2128 0.000275847 1.14736e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.21
## p 0.75
## se 0.73
## sp 0.22
## ppv 0.73
## npv 0.22mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -17.49 -17.17 1.33 1.11 -20.03 -16.04 1.01 745 996
## p 0.51 0.52 0.29 0.37 0.07 0.95 1.01 662 632
## se 0.69 0.71 0.13 0.13 0.47 0.88 1.00 2735 1625
## sp 0.27 0.25 0.12 0.12 0.09 0.50 1.01 1882 1217
## ppv 0.50 0.51 0.29 0.38 0.06 0.95 1.01 675 653
## npv 0.46 0.45 0.30 0.39 0.03 0.95 1.01 748 786
2x2 cross table
Effect occur y=1 with probabilty p01, p11 from each cause x{0,1}
event frequncy nxy of effect y{0,1} by cause x{0,1}
n01~B(n0.,p0)
n11~B(n1.,p1)
n01=n0p0, n00=n0(1-p0)
n11=n1p1, n10=n1(1-p1)
p00=n00/n=n0(1-p0)/n, p01=n01/n=n0p0/n
p10=n10/n=n1(1-p1)/n, p11=n11/n=n1p1/n
Cramer'V (chi2/n/(min(row,column)-1))^.5
in 2x2
crv =(n11*n00-n10*n01)/(n0.*n1.*n.0*n.1)^.5
=(n0(1-p0)n1p1-n0p0n1(1-p1))/(n0n1(n0(1-p0)+n1(1-p1))(n0p0n1p1))^.5
kappa coefficient k=(po-pe)/(1-pe)
po: Observed agreement (proportion of times both raters agreed)
pe: Expected agreement under independence
po=p00+p11
=(n0(1-p0)+n1p1)/n
pe=(p00+p01)(p00+p10)(p10+p11)(p01+p11)
=n0/n*(n0(1-p0)+n1(1-p1))/n*(n0p0+n1p1)/n*n1/nex16-1.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
}ex16-2.stan
data {
int n;
int n0;
int n01;
int n1;
int n11;
}
parameters {
real<lower=0,upper=1> p0;
real<lower=0,upper=1> p1;
}
model {
n01~binomial(n0,p0);
n11~binomial(n1,p1);
}
generated quantities {
real RR;
RR=p1/p0;
real OR;
OR=(p1/(1-p1))/(p0/(1-p0));
real crv;
crv=(n0*(1-p0)*n1*p1-n0*p0*n1*(1-p1))/(n0*n1*(n0*(1-p0)+n1*(1-p1))*(n0*p0+n1*p1))^.5;
real k;
real po;
po=(n0*(1-p0)+n1*p1)/n;
real pe;
pe=n0/n*(n0*(1-p0)+n1*(1-p1))/n*(n0*p0+n1*p1)/n*n1/n;
k=(po-pe)/(1-pe);
}n0=30
n01=rbinom(1,n0,0.3)
n1=30
n11=rbinom(1,n1,0.6)
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -44.3649
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -38.8529 0.000121331 1.41282e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -38.85
## p0 0.30
## p1 0.57
## RR 1.89
## OR 3.05mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -42.82 -42.51 0.97 0.74 -44.73 -41.86 1.01 981 1123
## p0 0.31 0.31 0.08 0.08 0.18 0.45 1.01 2044 1378
## p1 0.56 0.56 0.08 0.09 0.42 0.70 1.00 2172 1722
## RR 1.95 1.83 0.67 0.54 1.11 3.24 1.00 1972 1454
## OR 3.40 2.93 1.94 1.48 1.21 7.12 1.00 2062 1553
data=list(n=n0+n1,n0=n0,n01=n01,n1=n1,n11=n11)
mdl=cmdstan_model('./ex16-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -47.0484
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -38.8529 0.00190612 0.00050074 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -38.85
## p0 0.30
## p1 0.57
## RR 1.89
## OR 3.05
## crv 0.27
## k 0.61
## po 0.63
## pe 0.06mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -42.82 -42.51 0.97 0.74 -44.73 -41.86 1.01 981 1123
## p0 0.31 0.31 0.08 0.08 0.18 0.45 1.01 2044 1378
## p1 0.56 0.56 0.08 0.09 0.42 0.70 1.00 2172 1722
## RR 1.95 1.83 0.67 0.54 1.11 3.24 1.00 1972 1454
## OR 3.40 2.93 1.94 1.48 1.21 7.12 1.00 2062 1553
## crv 0.26 0.26 0.12 0.12 0.05 0.45 1.00 2082 1633
## k 0.60 0.60 0.06 0.06 0.49 0.70 1.00 2093 1659
## po 0.63 0.63 0.06 0.06 0.52 0.72 1.00 2092 1659
## pe 0.06 0.06 0.00 0.00 0.06 0.06 1.00 1809 1386
point of subjective equality PSE
PSE: 50% threshold for sensing the difference between two stimuli is equal
JND: Just noticeable difference, difference between 50% threshold and 75%
r~B(n,p) #reaction for stimuli
p=1/(1+exp(-(a+b*x)))
x=x1-x0, x0,x1 is strength of stimuli
PSE=-a/b
JND=(log(0.75/0.25)-a)/b-PSEex6-3-0.stan
mulit logistic regression
data{
int N;
int m;
vector[N] x;
array[N] int y;
}
parameters{
real b0;
real b1;
}
model{
y~binomial_logit(m,b0+b1*x);
}n=20
m=10
x=runif(n,-2,2)
y=rbinom(n,m,1/(1+exp(-(-2+3*x))))
glm(y/m~x,family=binomial('logit'))##
## Call: glm(formula = y/m ~ x, family = binomial("logit"))
##
## Coefficients:
## (Intercept) x
## -2.50 3.29
##
## Degrees of Freedom: 19 Total (i.e. Null); 18 Residual
## Null Deviance: 13.9
## Residual Deviance: 1.66 AIC: 13.5data=list(N=n,m=m,x=x,y=y)
mdl=cmdstan_model('./ex6-3-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -224.264
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -77.1961 0.00156079 2.52577e-06 1 1 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -77.20
## b0 -2.50
## b1 3.29mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -78.18 -77.87 1.00 0.70 -80.18 -77.25 1.00 748 767
## b0 -2.58 -2.56 0.45 0.44 -3.35 -1.88 1.01 526 489
## b1 3.38 3.37 0.47 0.47 2.62 4.17 1.01 535 552
b0=mcmc$draws('b0') |> as.vector()
b1=mcmc$draws('b1') |> as.vector()
pse=-b0/b1
quantile(pse,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.572 0.645 0.665 0.721 0.760 0.802 0.866 0.880 0.964jnd=(log(0.75/0.25)-b0)/b1-pse
quantile(jnd,probs=c(0.0,0.025,0.05,0.25,0.5,0.75,0.95,0.975,1),na.rm=T)## 0% 2.5% 5% 25% 50% 75% 95% 97.5% 100%
## 0.219 0.256 0.263 0.298 0.326 0.359 0.420 0.442 0.547x1=runif(length(b0),-2,2)
p=1/(1+exp(-(b0+b1*x1)))
pm=1/(1+exp(-(median(b0)+median(b1)*x1)))
qplot(x1,pm,col=I('darkgray'),ylab='p')+
geom_line(aes(x=x1,p=p),col='red')